function vtc = importvtcfromanalyze(imgs, bbox, res, imeth)
% importvtcfromanalyze  - import a VTC from Analzye files
%
% FORMAT:       vtc = importvtcfromanalyze(imgs [, bb [, res [, imeth]]])
%
% Input fields:
%
%       imgs        cell array with HDR filenames for BVQXfile
%       bb          optional 2x3 bounding box (default: MNI covering box)
%       res         optional resolution (default: 3)
%       imeth       interpolation, 'cubic', {'linear'}, 'nearest'
%
% Output fields:
%
%       vtc         created VTC object
%
% Note: this function requires the MEX file flexinterpn.

% Version:  v0.8a
% Build:    9102718
% Date:     Oct-27 2009, 6:14 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
   ~iscell(imgs) || ...
    numel(imgs) < 1 || ...
   ~ischar(imgs{1}) || ...
    isempty(imgs{1}) || ...
    exist(imgs{1}(:)', 'file') ~= 2
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad input argument.' ...
    );
end

% go on with loading images and sampling VTC data
imgs = imgs(:);
nimg = numel(imgs);
himg = cell(1, nimg);
vimg = zeros(20, nimg);
nvol = ones(1, nimg);
try
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', 'Converting Analyze to VTC...');
    BVQXprogress(pbar, 0, 'Checking HDRs...', 'visible', 0, 8 * nimg);
catch
    pbar = [];
end

% disable unwindstack for now...
uws = BVQXfile(0, 'unwindstack');
BVQXfile(0, 'unwindstack', false);
try
    for ic = 1:nimg
        if ~ischar(imgs{ic}) || ...
            isempty(imgs{ic})
            error('BAD_IMAGENAME');
        end
        imgs{ic} = strrep(strrep(imgs{ic}, '.img', '.hdr'), '.IMG', '.HDR');
        himg{ic} = BVQXfile(imgs{ic});
        if numel(himg{ic}) ~= 1 ||...
           ~isBVQXfile(himg{ic}, 'hdr')
            error('BAD_IMAGECONT');
        end
        icf = himg{ic}.CoordinateFrame;
        vimg(:, ic) = [icf.Trf(:); icf.Dimensions(:)];
        nvol(ic) = size(himg{ic}.VoxelData, 4);
        if ~isempty(pbar)
            BVQXprogress(pbar, ic);
        end
    end
    if any(any(diff(vimg, 1, 2)))
        warning( ...
            'BVQXtools:BadArgument', ...
            'Spatial orientation/dimensions mismatch between images.' ...
        );
    end
catch
    clearbvqxobjects(himg);
    BVQXfile(0, 'unwindstack', uws);
    if ~isempty(pbar)
        closebar(pbar);
    end
    error( ...
        'BVQXtools:BadArgument', ...
        'Error loading image %d.', ...
        ic ...
    );
end

% checking other arguments
sfn = himg{1}.FilenameOnDisk;
if nargin < 2 || ...
   ~isa(bbox, 'double') || ...
   ~isequal(size(bbox), [2, 3]) || ...
    any(isnan(bbox(:)) | bbox(:) < 0 | bbox(:) > 255)
    bbox = [];
else
    bbox = round(bbox);
end
if nargin < 3 || ...
   ~isa(res, 'double') || ...
    numel(res) ~= 1 || ...
   ~any((1:15) == res)
    res = 3;
end
if nargin < 4 || ...
   ~ischar(imeth) || ...
   ~any(strcmpi(imeth(:)', {'cubic', 'linear', 'nearest'}))
    imeth = 'linear';
else
    imeth = lower(imeth(:)');
end

% get global setting to figure out DataType/FileVersion
global bvqxfile_config;
dtype = bvqxfile_config.settings.DataTypes.VTC;    

% guess scaling
if dtype == 1
    if istransio(himg{1}.VoxelData)
        vdm = minmaxmean(himg{1}.VoxelData(:, :, :));
    else
        vdm = minmaxmean(himg{1}.VoxelData);
    end
    vdm = vdm(2);
    if vdm > 16384
        vdf = 16384 / double(vdm);
    else
        vdf = [];
    end
else
    vdf = 1;
end

% try to sample first vol (check)
try
    [vtd, obox] = himg{1}.SampleBox(struct('BBox', bbox, 'XYZRes', res), 1, imeth);
    if dtype == 1
        vtd = shiftdim(uint16(round(vtd)), -1);
    else
        vtd = shiftdim(single(vtd), -1);
    end
catch
    clearbvqxobjects(himg);
    BVQXfile(0, 'unwindstack', uws);
    rethrow(lasterror);
end

% create VTC data
vtd(sum(nvol), 1, 1, 1) = 0;
vts = size(vtd);

% sampling
if ~isempty(pbar)
    BVQXprogress(pbar, nimg, 'Sampling images...');
end
icc = 1;
for ic = 1:nimg
    for vc = 1:nvol(ic)
        try
            hy = shiftdim(himg{ic}.SampleBox( ...
                struct('BBox', bbox, 'XYZRes', res), vc, imeth), -1);
            if ~isempty(vdf)
                hy = vdf * hy;
            end
            if ~isempty(pbar)
                BVQXprogress(pbar, nimg + 7 * ((ic - 1) + vc / nvol(ic)));
            end
            vtd(icc, :, :, :) = hy;
            icc = icc + 1;
        catch
            clearbvqxobjects(himg);
            BVQXfile(0, 'unwindstack', uws);
            if ~isempty(pbar)
                closebar(pbar);
            end
            error( ...
                'BVQXtools:InternalError', ...
                'Error sampling data of volume %d.', ...
                ic ...
            );
        end
    end
    himg{ic}.ClearObject;
    himg{ic} = [];
end

% create VTC
vtc = bless(BVQXfile('new:vtc'), 1);
vtc.DataType = dtype;
if dtype > 1
    vtc.FileVersion = 3;
end
vtc.NameOfSourceFMR = sfn;
vtc.Resolution = res;
vtc.XStart = obox.BBox(1, 1);
vtc.XEnd = obox.BBox(1, 1) + res * vts(2);
vtc.YStart = obox.BBox(1, 2);
vtc.YEnd = obox.BBox(1, 2) + res * vts(3);
vtc.ZStart = obox.BBox(1, 3);
vtc.ZEnd = obox.BBox(1, 3) + res * vts(4);
vtc.VTCData = vtd;
vtc.NrOfVolumes = size(vtc.VTCData, 1);

% close bar
if ~isempty(pbar)
    closebar(pbar);
end

% reset unwindstack flag
BVQXfile(0, 'unwindstack', uws);
